#ifndef HMS_CUBIC_ROOT_HPP
#define HMS_CUBIC_ROOT_HPP

#include "globals.hpp"
#include "unifiedMathFunctions.hpp"
#include "typeAliases.hpp"

#include <cmath>

namespace hms::cubic
{


/**
 * @brief Returns a coeff. of a depressed cubic function
 * 
 * @param a coeff of a cubic function
 * @param b coeff of a cubic function
 * @param c coeff of a cubic function
 */
template<typename A, typename B, typename C>
auto p( const A& a, const B& b, const C& c ) -> decltype(auto) {
	return ( 3 * a * c - square(b) ) / ( 3 * square(a) );
}


/**
 * @brief Returns a coeff. of a depressed cubic function
 * 
 * @param a coeff of a cubic function
 * @param b coeff of a cubic function
 * @param c coeff of a cubic function
 * @param d coeff of a cubic function
 */
template<typename A, typename B, typename C, typename D>
auto q( const A& a, const B& b, const C& c, const D& d ) -> decltype(auto) {
	return ( 2 * cube(b) - 9*a*b*c + 27*square(a) * d ) / ( 27 * cube(a) );
}

/**
 * @brief Returns the discriminant D of a depressed cubic.
 * D < 0 -> three real roots
 * D > 0 -> one real root
 * D = 0, p != 0 -> one single, one double root
 * D = 0, p = 0 -> triple root at 0
 * 
 * @tparam P 
 * @tparam Q 
 * @param p coeff of depressed cubic
 * @param q coeff of depressed cubic
 * @return decltype(auto) 
 */
template<typename P, typename Q>
auto discriminant( const P& p, const Q& q ) -> decltype(auto) {
	return 4 * cube(p) + 27 * square(q);
}


/**
 * @brief Converts the root of a depressed cubic into that of the original 
 * cubic function.
 * 
 * @tparam T 
 * @tparam A 
 * @tparam B 
 * @param depressedRoot root of the depressed cubic to be converted
 * @param a coeff of a cubic function
 * @param b coeff of a cubic function
 * @return decltype(auto) 
 */
template<typename T, typename A, typename B>
auto root( const T& depressedRoot, const A& a, const B& b) -> decltype(auto) {
	return depressedRoot - (b / (3*a));
}

/**
 * @brief In the case of a positive cubic::discriminant,
 * the cubic equation has a single real root,
 * and the functions in this namespace apply.
 * 
 */
namespace singleReal
{

/**
 * @brief Helper variable for the single real root case:
 * s = sqrt( (q/2)^2 + (p/3)^3 )
 * Needed for the two summands of Cardano's original formulation for the root r:
 * r = C1 + C2
 * C1 = (-q/2 + s)^(1/3)
 * C2 = (-q/2 - s)^(1/3)
 * 
 * @tparam P 
 * @tparam Q 
 * @param p coeff of depressed cubic
 * @param q coeff of depressed cubic
 * @return auto 
 */
template<typename P, typename Q>
auto sqrtTerm( const P& p, const Q& q ) -> decltype(auto) {
	using Eigen::sqrt;
	using std::sqrt;
	return sqrt( 0.25*square(q) + cube(p) / 27. );
}

template<bool negative, typename Q, typename S>
auto cubeRoot( const Q& q, const S& sqrtTerm ) -> decltype(auto) {
	static constexpr int sign { negative ? -1 : 1 };
	return pow( -q/2 + sign * sqrtTerm, 1./3 );
}

template<typename Q, typename S>
auto depressedRoot( const Q& q, const S& sqrtTerm ) -> decltype(auto) {
	return cubeRoot<0>(q, sqrtTerm) + cubeRoot<1>(q, sqrtTerm);
}

} // namespace singleReal

/**
 * @brief In the case of a negative cubic::discriminant (or zero, with p <> 0),
 * the cubic equation has three real roots.
 * Then the functions in this namespace apply.
 * 
 */
namespace threeReal
{

/**
 * @brief Returns the scale factor 2 * sqrt(-p/3)
 * for the trig. computation of cubic roots
 * 
 * @tparam P 
 * @param p coeff of depressed cubic
 * @return decltype(auto) 
 */
template<typename P>
auto scale( const P& p ) -> decltype(auto) {
	using Eigen::sqrt;
	using std::sqrt;
	return 2 * sqrt( - p / 3 );
}

/**
 * @brief Returns the acos-Term 
 * 1/3 * acos( 1.5 q/p * sqrt(-3/p) )
 * for the trig. computation of cubic roots
 * 
 * @param p coeff of depressed cubic
 * @param q coeff of depressed cubic
 * @return decltype(auto) 
 */
template<typename P, typename Q>
auto acosTerm( const P& p, const Q& q ) -> decltype(auto) {
	using Eigen::acos;
	using Eigen::sqrt;
	using std::acos;
	using std::sqrt;
	return 1./3 * acos( 1.5 * q / p * sqrt(-3 / p) );
}

/**
 * @brief Returns one of three depressed cubic roots, depending on k = {0,1,2}.
 * The equation is:
 * 
 * scale * cos(acosTerm - k * 2 pi / 3)
 * 
 * @tparam S 
 * @tparam A 
 * @tparam K 
 * @param scale
 * @param acosTerm
 * @param k 
 * @return decltype(auto) 
 */
template<typename S, typename A, typename K>
auto depressedRoot( const S& scale, const A& acosTerm, const K& k ) -> decltype(auto) {
	using Eigen::cos;
	using std::cos;
	return scale * cos( acosTerm - k * 2 * math::pi / 3 );
}

} // namespace threeReal

/**
 * @brief Cubic equation which is normalised and has no linear term, i.e.
 * ax^3 + bx^2 + cx + d = 0,
 * where a = 1 and c = 0, therefore
 * x^3 + bx^2 + d = 0
 * 
 * If instead you have the form
 * x^3 + cx + d,
 * then it is equivalent to the depressed form
 * t^3 + pt + q,
 * i.e. c = p and d = q.
 * Therefore, you can check the discriminant with discriminant(p,q) and then
 * use threeReal::scale and threeReal::acosTerm to compute the threeReal::depressedRoot,
 * which is the actual root in this case.
 */
namespace norm_noLin
{

/**
 * @brief Returns coeff p of a depressed cubic,
 * in the simplified case of a normalised cubic function without a linear term,
 * i.e. a = 1, c = 0.
 * 
 * @see hms::cubic::p
 */
template<typename B>
auto p( const B& b ) -> decltype(auto) {
	return - square(b) / 3;
}

/**
 * @brief Returns coeff q of a depressed cubic,
 * in the simplified case of a normalised cubic function without a linear term,
 * i.e. a = 1, c = 0.
 * 
 * @see hms::cubic::norm_noLin::p
 * @see hms::cubic::q
 */
template<typename B, typename D>
auto q( const B& b, const D& d ) -> decltype(auto) {
	return 2 * cube(b) / 27 + d;
}

template<typename T, typename B>
auto root( const T& depressedRoot, const B& b) -> decltype(auto) {
	return depressedRoot - (b / 3);
}

/**
 * @brief Number of independent helper variables per root
 * 
 */
constexpr Index bufCols {5};

/**
 * @brief Computation of cubic roots in the case of three real roots.
 * Takes parameters b and d of a cubic equation of the form
 * x^3 + bx^2 + d = 0,
 * as well as k = {0,1,2}, which decides which root to choose,
 * and returns the root x.
 * The discriminant is computed, to ascertain that it is negative or zero,
 * in which case three real roots exist. Otherwise, an exception is thrown.
 * For k, see the plot below, which assumes three real roots:
 *      
 *      ^ y
 *      |         /
 *     _|        /
 * ___/_|\______/_______\ 
 *  / ^ | ^\__/ k=0     / x
 * /  | | |k=1
 *  k=2
 * 
 * @tparam B 
 * @tparam D 
 * @tparam K 
 * @tparam Buf 
 * @param bXpr Parameter b, preferrably as Eigen expression
 * @param dXpr Parameter d, preferrably as Eigen expression
 * @param k root selection parameter
 * @param calcBuf Buffer for helper variables.
 * Must be at least the size of b times bufCols
 * @return decltype(auto) 
 */
template<typename B, typename D, typename K, typename Buf>
auto root(
	const B& bXpr, const D& dXpr,
	const K& k, 
	Buf& calcBuf
) -> decltype(auto);


/**
 * @brief Computation of cubic roots in the case of a single real root.
 * Takes parameters b and d of a cubic equation of the form
 * x^3 + bx^2 + d = 0,
 * and returns the root x.
 * The discriminant is computed, to ascertain that it is positive,
 * in which case a single real root exists. Otherwise, an exception is thrown.
 * 
 * @tparam B 
 * @tparam D 
 * @tparam K 
 * @tparam Buf 
 * @param bXpr Parameter b, preferrably as Eigen expression
 * @param dXpr Parameter d, preferrably as Eigen expression
 * @param calcBuf Buffer for helper variables.
 * Must be at least the size of b times bufCols
 * @return decltype(auto) 
 */
template<typename B, typename D, typename Buf>
auto root(
	const B& bXpr, const D& dXpr,
	Buf& calcBuf
) -> decltype(auto);


} // namespace norm_noLin


} // namespace hms::cubic

#include "cubicRoot.tpp"

#endif